Single-lead ECG based autonomic nervous system assessment for meditation monitoring

We propose a single-lead ECG-based heart rate variability (HRV) analysis algorithm to quantify autonomic nervous system activity during meditation. Respiratory sinus arrhythmia (RSA) induced by breathing is a dominant component of HRV, but its frequency depends on an individual’s breathing speed. To address this RSA issue, we designed a novel HRV tachogram decomposition algorithm and new HRV indices. The proposed method was validated by using a simulation, and applied to our experimental (mindfulness meditation) data and the WESAD open-source data. During meditation, our proposed HRV indices related to vagal and sympathetic tones were significantly increased (p < 0.000005) and decreased (p < 0.000005), respectively. These results were consistent with self-reports and experimental protocols, and identified parasympathetic activation and sympathetic inhibition during meditation. In conclusion, the proposed method successfully assessed autonomic nervous system activity during meditation when respiration influences disrupted classical HRV. The proposed method can be considered a reliable approach to quantify autonomic nervous system activity.

ECG based cardiopulmonary coupling analysis and HRV analysis have been used to assess psychiatric conditions. Cardiopulmonary coupling analysis measures the synchronization between the heartbeat interval and respiration by using the spectral coherence, cross entropy and phase locking value (PLV) [20][21][22] . This synchronization phenomenon is called respiratory sinus arrhythmia (RSA), which is the modulation of efferent parasympathetic nerves caused by baroreflexes (arterial baroreflex, lung stretch reflex, and Bainbridge reflex) during the inspiration phase [23][24][25][26] . RSA is affected by the tidal volume and respiratory rate (RR) [27][28][29] , and its quantity reflects sleep quality and apnea 20,22 . To capture other neurophysiological changes, HRV is widely used, and its frequency analysis reflects the balance between vagal and sympathetic tones. The former shows a fast response, and the latter dominates the low-frequency band of HRV and seems to be independent of the baroreflex 30 .
Meditation (or relaxation) training leads to parasympathetic activation and sympathetic inhibition, and it is supported by several modalities, such as muscle sympathetic nerve activity 31 , functional MRI 2 , and electroencephalogram 32 . Some HRV studies on meditation also reported an increase in high-frequency power [32][33][34][35] , which indicates vagal tone, but other studies showed the opposite results [36][37][38][39][40] . These contradictory results seem to be due to the influence of RSA. When the RR < 0.15 Hz, the RSA component moves from the high-to low-frequency region, and then the low-and high-frequency powers can be overestimated and underestimated, respectively 28 . Hence, the fluctuation of an individual's RR led to significant variation in HRV analysis 41 , and there was no report on the identification of cardiac sympathetic inhibition during meditation.
To separate the RSA effect and other autonomic nervous system activities in HRV analysis, various decomposition algorithms have been utilized, such as independent component analysis 42 , adaptive noise cancellers 43 , autoregressive moving average with exogenous input model (ARMAx) 44 , and orthogonal subspace projection (OSP) 45,46 . The adaptive noise canceller, ARMAx, and OSP share the same algorithmic structure, with OSP exhibiting the best performance 45 . Even though these algorithms significantly enhance the reliability of HRV, they assume a linear relationship between the RSA component and a respiration signal, while in reality the relationship is nonlinear 47 . As an alternative, a notch filter with a Gaussian bell shape (Gauss) was employed, but it assumed that the instantaneous RR is constant 48 . In this study, to overcome the nonlinearity and nonstationarity of RSA, we proposed a novel HRV tachogram decomposition algorithm, and new HRV indices to quantify vagal and sympathetic tones during meditation (see Fig. 2).
The proposed method is described in the "Materials and methods" section. We compared the performance of decomposition algorithms using simulation, and observed the changes in the vagal and sympathetic tones during meditation through real ECG data in "Results" section. In "Discussion" section, the results are discussed.

Materials and methods
Since HRV reflects the balance of the autonomic nervous system, HRV analysis is frequently used to assess stress and psychiatric conditions. The first step of HRV analysis is QRS detection on the ECG. In this study, we bandpass filtered an ECG with 0.5-50 Hz by a 3rd-order Butterworth, and we employed the Pan-Tompkins algorithm to detect R-peaks of the ECG 49 . The interval between successive R-peaks (RRI) is accumulated over a predefined   www.nature.com/scientificreports/ time window. In this study, we adopted a 5 min time window for short-term HRV analysis. To detect ectopic beats and incorrect peaks, we calculated the k-th derivative of the instantaneous heart rate (HR) 50 as follows: where R k is the k-th R-peak. When an HR' k was greater than a predefined threshold, R k was excluded from the HRV analysis, and the corrected R-peaks were double-checked by human experts. Because the RRI is sampled at irregular time points (heartbeats), the RRI must be converted into an evenly sampled time series for frequency analysis. A uniformly sampled signal of the RRI, called a tachogram, was reconstructed using cubic spline interpolation, and we resampled the tachogram to 4 Hz, as in previous studies [44][45][46] . For precise HRV analysis, it is necessary to distinguish the cardiac vagal tone modulated by respiration and other autonomic nervous system activities. Hence, we consider the original tachogram x RAW as follows: where x RSA and x R-f are the RSA and RSA-free tachograms, respectively.
x R-f is obtained by removing x RSA from x RAW . The tachogram decomposition algorithm will be described in the next subsection. In this study, since we focus on the tachogram decomposition algorithm, we adopt only frequency domain HRV indices.
ECG-derived respiration. To measure the respiration signal, various sensors have been utilized, such as a pressure transducer and a thermistor. These sensors provide accurate information but are both expensive and inconvenient to wear. As an alternative, several algorithms to estimate the respiration signal from biomedical signals, such as photoplethysmogram and ECG, have been proposed [51][52][53][54] . In this study, we employed an ECGderived respiration (EDR) algorithm that captures the change in electrical impedance and the cardiac vector according to chest expansion. Here, we briefly review 4 frequently used EDR algorithms: EDR A (R-wave amplitude), EDR QR (Q-R slope), EDR RS (R-S slope) and EDR SR (slope range). Each EDR algorithm uses a respirationrelated sample r(i) at the i-th heartbeat as follows: • EDR A : r A (i) = R-peak amplitude • EDR QR : r QR (i) = steepest ascent slope on Q-R wave • EDR RS : r RS (i) = steepest descent slope on R-S wave • EDR SR : r SR (i) = r QR (i) − r RS (i) By interpolating the respiration-related samples, the EDR signal is generated. The overall procedure is explained in detail in 51 .

Tachogram decomposition. Previous algorithms.
To decompose raw tachogram x RAW into x RSA and x R-f , various algorithms have been proposed, such as Gauss, ARMAx and OSP. Both the ARMAx and OSP algorithms assume x RSA as follows: where RESP is a respiratory basis matrix. The trajectory matrix 46 or wavelet coefficients 45 of EDR have been utilized as the RESP matrix. In this study, the trajectory matrix was employed for the RESP. We assigned 3 s to the lag size of the trajectory matrix of the ARMAx, while the lag size of OSP was set to the minimum value between the minimum description length and the Akaike information criterion, with a maximum lag of 10 s. The RSA-free tachogram is estimated as follows: The ARMAx and the OSP are described in [44][45][46] in detail. As another tachogram decomposition algorithm, a notch filter with Gaussian bell shape (Gauss) was proposed. Its magnitude-squared frequency response is: where θ represents the notch frequency and corresponds to the average RR. α and β are fitting parameters and it can be computed from the spectrum of x RAW by using a least square method. The filtered signal will be x R−f , and x RSA is acquired by subtracting x R−f from x RAW . The overall procedure of Gauss is described in 48 .
Proposed algorithm. Based on the synchronization between the frequency of x RSA and the instantaneous RR, we designed a novel tachogram decomposition algorithm. Because it has a zero-phase response, we call it the zero-phase line enhancer (ZLE). The ZLE can be interpreted as an adaptive filter without phase distortion 55 . As shown in Fig. 3, the ZLE has two parts: (a) instantaneous RR estimation and (b) time-variant forward-backward IIR notch filtering. To estimate the instantaneous RR, we employed a smoothed pseudo Wigner-Ville distribution of the EDR signal 56,57 as follows: where r is the pole-zero contraction factor, which determines the stop band. When r is closer to 1, the frequency response of the notch filter becomes narrower. We set r to 0.93, and then the stop band of H(z, θ[n]) is 0.047 Hz. As shown in Fig. 3, the time-variant forward-backward IIR notch filtering performs four sequential steps: 1-forward IIR notch filtering, 2-time reversal, 3-backward IIR notch filtering, and 4-time reversal. The zero-phase response can be explained in the z-domain as follows: www.nature.com/scientificreports/ where F(z) and B(z) mean forward-and backward-filtered signals, respectively. P(z) represents the time-reversed F(z). The filter output is the estimated RSA-free tachogram x R−f , and the RSA tachogram is estimated as follows: RSA indices. RSA has been utilized as a measure of cardiac vagal function. RSA can be quantified by cardiopulmonary coupling analysis. In this study, we utilized the PLV 22,39,59 . The PLV is broadly used to measure a phase interaction between two signals and is defined as follows: where φ 1 (n) and φ 2 (n) represent the phases of the respiration signal (or EDR) and tachogram x RAW , respectively. The phases are calculated by a Hilbert transform. As another measure for RSA, the extent of RSA (RSA E ) 27,39,60 was calculated by the median value of RRI variation (longest RRI-shortest RRI) within each breathing cycle as follows: where RRI k represents the subset of RRIs within the k-th breathing cycle. max and min are maximum and minimum functions, respectively. K indicates the number of breaths taken in the time window (5 min). Because RSA E is based on the RRI, the unit of RSA E is ms.

HRV indices.
In HRV analysis, the power spectrum of x RAW is segmented into low-frequency (0.04-0.15 Hz) and high-frequency (0.15-0.4 Hz) bands, and each band power is called LF (low frequency) and HF (high frequency) 61 . HF is regarded as an index for vagal tone, and LF reflects both sympathetic and vagal tones. To quantify sympathetic tone, LF/HF is generally used. Since these HRV indices were insufficient to cover the aforementioned RSA issue, it is necessary to simultaneously use both power spectra of x R−f and x RSA ; Subscripts R-f and RSA were used to distinguish each band power ( x R−f : LF R-f and HF R-f , x RSA : LF RSA and HF RSA ). To assess sympathovagal balance, Varon et al. designed an HRV index SB U based on the separated band powers 46 as follows: Since x RSA and HF R-f purely indicate vagal tone but LF R-f is related to both vagal and sympathetic tones, we proposed the new HRV indices related to vagal (Eq. 17) and sympathetic (Eq. 18) tones as follows: We employed Welch's method to compute the power spectra.

Simulation.
To compare the tachogram decomposition algorithms, we designed a simulated tachogram x RAW . First, to cover various situations, we considered three respiration types as follows: where ψ (n) , the derivative of ψ(n) , means the simulated RR. ①, ②, and ③ represent the frequencies of a sine wave (0.15 Hz), linear chirp (0.04-0.4 Hz), and sinusoidal frequency modulated signal (0.04-0.4 Hz). Since not only the frequency but also the amplitude of the tachogram is modulated by RR 27-29 , we devised an amplitude modulation function for the simulated RSA component x RSA . Specifically, because the amplitude of the www.nature.com/scientificreports/ tachogram with respect to RR in 29 was skewed to the right, we utilized the gamma distribution function as the amplitude modulation function m() as follows: where Г(α) is the gamma function. α and β are shape and rate parameters (α = 2, β = 15). Note that when ψ (n) is constant the amplitude modulation function m() becomes constant. On the other hand, as the simulated RSAfree tachogram x R-f , we employed white Gaussian noise as follows: where σ 2 is the variance of the white Gaussian noise. We adjusted σ 2 so that the power ratio of x RSA and x R-f was equal (0 dB), and then the simulated tachogram x RAW was given by their summation.
For rigorous comparison, we performed Monte Carlo simulation with 1000 independent simulations.
Database. To investigate the changes in HRV indices according to the tachogram decomposition algorithms during meditation, we utilized two real ECG databases: mindfulness meditation data and WESAD data.
Mindfulness meditation data. We collected a single-lead ECG and a respiration signal from 8 male and 8 female subjects (age = 28.5 ± 4.2 years) during mindfulness meditation. All subjects were novice meditators who had engaged in meditation within 1-2 times. All subjects rested for 10 min and then performed mindfulness meditation for 40 min under the guidance of a meditation teacher. During mindfulness meditation, subjects were guided to maintain awareness of their breathing and sensations without judgment. Breathing speed was not restricted. All subjects closed their eyes and remained in a sitting position to minimize motion artifacts in the ECG. They completed a visual analog scale (VAS) for stress, tension, and concentration before and after mindfulness meditation. Each ECG was measured by a wearable ECG patch (EP200, Life Science Technology ® , Seoul, Korea), and the sampling rate was 250 Hz. The reference respiration signal was recorded by a respiratory transducer (RSPEC-R, Biopac ® , California, USA) with a sampling rate of 250 Hz. We used the reference respiration signal to evaluate EDR and measure vital signs and RSA indices, and we employed EDR for the tachogram decomposition algorithms. Informed consent was obtained from all subjects, and the institutional review board of the Korea Institute of Science and Technology (KIST) approved all procedures in this study (KIST-2019-015, July. 19, 2019). Informed consent was obtained from all participants prior to the experiment, and all experiments were conducted in strict accordance with KIST ethics guidelines and the declaration of Helsinki.

WESAD data.
We used open-source data (WESAD, Wearable Stress and Affect Detection) to verify the proposed method in a public domain. This dataset contains the ECGs of 15 subjects (age = 27.5 ± 2.4 years), and its sampling rate is 700 Hz. The experimental protocol consisted of four different affective states (neutral, stress, amusement, meditation), and each emotion was mapped onto a two-dimensional affective space (arousal and valence) by using self-reports (SAM, Self-Assessment Manikins). The WESAD data are described in detail in 62 .
Statistical evaluation. Since the simulation was performed with white Gaussian noise, ANOVA and t tests were used to compare the tachogram decomposition algorithms. In the case of real data (mindfulness meditation and WESAD data), nonparametric statistical tests were utilized because the number of subjects was not large. Specifically, the Wilcoxon signed rank test was employed when the two groups had the same sample size, and the Mann-Whitney U test was used otherwise. We defined p < 0.000005 as statistically significant. To depict the distribution, we utilized a box plot where the top and bottom boxes represent the 75th and 25th percentiles, and the center, top, and bottom lines represent the 50th, 90th, and 10th percentiles, respectively.

Results
Simulation. For each simulation type (RR: constant, linear function, and sinusoidal function), we calculated the Pearson correlation coefficient and root mean square error (RMSE) between x R-f and x R-f of the tachogram decomposition algorithms (Gauss, ARMAx, OSP, and ZLE). Figure 4 depicts the box plots of correlation coefficients and RMSEs ((a): a sine wave, (b): linear chirp, and (c): sinusoidal frequency modulated signal). For each simulation type, we conducted one-way ANOVA with Bonferroni correction to compare four algorithms (Gauss, ARMAx, OSP, and ZLE), and significant differences (p < 0.000005) were observed in both the Pearson correlation coefficients and RMSEs for all simulation types. For pairwise comparisons, post hoc tests were performed by using paired t tests, and all algorithm comparisons showed significant differences (p < 0.000005). The OSP showed the lowest RMSE and highest correlation coefficient when the simulated RR was constant (see Fig. 4a), but the ZLE had the lowest RMSE and highest correlation coefficient for other RR types (see Fig. 4b,c). It is notable that ZLE was robust to the fluctuation of instantaneous RR but contrastively Gauss did not work well in time-varying RR.
(21) x RSA (n) = m(ψ(n)) · cos(2πψ(n)) EDR SR ) by using the reference respiration signal of the mindfulness meditation data. Because there was a phase difference between EDR and reference respiration signal, we utilized two measures that were not significantly affected by the phase difference. One is the PLV between each EDR and the reference respiration signal, and the other is the mean absolute error (MAE) between RRs calculated by each EDR and the reference respiration signal. The PLV and MAE reflect the similarity of phases and the accuracy of RR estimation, respectively. As shown in Table 1, EDR SR showed the best performance as in previous study 51 , so we employed EDR SR in the following analysis.
Mindfulness meditation data. Before applying the tachogram decomposition algorithms, we observed the effect of mindfulness meditation through various measures: VAS for stress, RSA indices, and vital signs. First, we confirmed the stress reduction after meditation through the VAS for stress (5.35 ± 1.77 → 2.47 ± 1.01). Second, we calculated the vital signs and RSA indices with 5 min intervals, and their box plots are depicted in To assess the trend of the HRV indices, we divided each record (50 min) into 1 rest phase and 4 meditation phases of equal length (10 min) and then compared the rest and each meditation phase by using a Wilcoxon's two-sample signed rank test (see Table 2). Table 2 depicts the mean and standard deviation of HRV indices. HF and IDX PNA are the HRV indices for vagal tone, and other indices (LF/HF, SB U and IDX SNA ) indicate sympathetic tone. No statistically significant differences were observed in HFs, but all IDX PNA s of ZLE were significantly increased during meditation (p < 0.000005). Unlike ZLE, IDX PNA s of previous tachogram decomposition algorithms did not significantly decrease at 20-30 min. For the HRV indices related to sympathetic tone, the LF/HF of the raw tachogram (x RAW ) significantly increased at 30-50 min, which seems to be caused by RSA, but all SB U s and IDX SNA s of ZLE were significantly diminished   Fig. 7; Table 2). Although, after 40 min, all SB U s and IDX SNA s of all algorithms except for Gauss were significantly decreased, but no statistically significant differences were observed in most SB U s and IDX SNA s of previous tachogram decomposition algorithms (Gauss, ARMAx, and OSP).
WESAD data. The WESAD data contains ECGs with 4 affective states: neutral, stress, amusement, and meditation. The SAM self-reports showed that the lowest arousal scores (2.3 ± 1.4) occurred in the meditation state. As shown in Fig. 8, the RSA indices (RSA E and PLV) and vital signs (HR and RR) had the highest and lowest values in the meditation state, respectively. For HR, RR, RSA E and PLV, we performed statistical tests between the neutral and meditation states by using Mann-Whitney U tests. Statistically significant differences (p < 0.000005) were observed in all comparisons except HR. This is similar to the results of the previous subsection. Figure Table 3. No statistically significant differences were observed in HFs, but IDX PNA of ARMAx and ZLE significantly increased during meditation, in parallel  www.nature.com/scientificreports/ with the RSA indices. Similar to the results of the previous subsection, SB U and IDX SNA of the ZLE significantly (p < 0.000005) decreased during meditation, but no statistically significant differences were observed in SB U and IDX SNA of other decomposition algorithms (Gauss, ARMAx and OSP) as shown in Fig. 9 and Table 3.

Discussion
The tachogram of HRV consists of various oscillations caused by complicated interactions between the heart and brain. In particular, the oscillation caused by RSA often forms most of the raw tachogram x RAW . For the precise quantification of autonomic nervous system activity, it is necessary to decompose the raw tachogram into RSA and RSA-free tachograms. In this study, we proposed the novel tachogram decomposition algorithm (ZLE) and new HRV indices (IDX PNA and IDX SNA ). We evaluated the performance of the proposed method through simulations, mindfulness meditation data, and WESAD data. RSA, which is a synchrony between the heartbeat interval and breathing, reflects cardiac vagal control. To quantify RSA, we employed the PLV and RSA E . Although the RSA E and PLV capture different features (phase interaction and RRI variation), these indices increased during meditation (see Figs. 6,9) as in previous study 39 , which indicates parasympathetic activation. In HRV analysis, all IDX PNA s of ZLE showed significant increases   Tables 2, 3), and IDX PNA appeared to be better than HF in identifying parasympathetic activation. Sympathetic inhibition during meditation is supported not only by literature but also by experimental analysis. Breathing meditation and slow breathing have been reported to improve pulmonary gas exchange performance 63 and lower blood pressure 64 , galvanic skin response 65 and muscle sympathetic nerve activity [66][67][68] . Reductions in vital signs (HR and RR) and self-report scores (stress and arousal) were also observed during meditation in our experiments, as shown in Figs. 6 and 9. Despite much evidence for sympathetic inhibition, the LF/HF of the classical HRV was noticeably increased during meditation [35][36][37][38][39][40] (see Tables 2, 3); it seems to be due to the influence of RSA. To cancel the interference of RSA, C. Varon et al. suggested SB U 46 and we proposed IDX SNA . Although SB U (LF R-f /(LF RSA + HF RSA )) is similar to IDX SNA (LF R-f /(LF RSA + HF RSA + HF R-f )), IDX SNA includes RSA-free vagal tone information (HF R-f ). Since RSA is insufficient to represent a whole cardiac vagal tone 69 , the proposed IDX SNA can be considered a more reasonable index to quantify sympathetic tone. In real meditation data, both SB U and IDX SNA of ZLE successfully identified sympathetic inhibition, and it may be the first report on the identification of cardiac sympathetic inhibition during meditation. As a result, both SB U and IDX SNA appear to be more reliable than classical LF/HF to identify sympathetic inhibition.
ARMAx and OSP share the same scheme, which finds a weight vector (w = (RESP T •RESP) −1 RESP•X) and subtracts x RSA from the original tachogram x RAW 45 . When the RESP is the trajectory matrix of EDR, the ARMAx and OSP can be interpreted as adaptive noise cancellers with a least square solution; then, the EDR and the RSA component (x RSA ) can be regarded as a noise reference signal and true noise, respectively. The algorithm performance is determined by a linear correlation between x RSA and the EDR 70 , but their true relationship is nonlinear 47 . Specifically, both the amplitude and frequency of the original tachogram are modulated by the RR    [27][28][29] , and such modulations cannot be expressed by a linear transform of EDR. As an alternative, Gauss was suggested, which assumes paced respiration 48 , but instantaneous RR generally changes with time. In the simulation, previous decomposition algorithms (Gauss, ARMAx, and OSP) performed well when RR was constant, but their performance was insufficient in time-varying RR. However, ZLE was robust to the fluctuations of simulated RR (see Fig. 4). Note that abrupt changes in RR often occur during breathing meditation. In fact, the HRV indices (IDX PNA , SB U , and IDX SNA ) of ZLE certainly identified both parasympathetic activation and sympathetic inhibition during meditation when no statistically significant differences were observed in some IDX PNA s and most IDX SNA of both ARMAx and OSP (see Tables 2, 3). Therefore, the ZLE seems to be the most reliable tachogram decomposition algorithm when RR dynamically changes. Despite the promising performance of the proposed ZLE and HRV indices, they have limitations. First, if some oscillation of the RSA-free tachogram (x R-f ) overlaps with the stop band of the notch filter in the frequency domain, then the corresponding oscillation will be suppressed. This filtering loss can be observed through the correlation coefficient between x RAW and x R-f when the influence of RSA was small (e.g., resting state). As shown in Fig. 6, the loss of ZLE seems to be insignificant. Second, the proposed method did not identify sympathetic activation. In the stress state, sympathetic nerves are activated, but none of the SB U s or IDX SNA s show a significant increase. In future works, we plan to advance the proposed method for stress assessment.

Conclusions
To precisely assess autonomic nervous system activity, we proposed a novel tachogram decomposition algorithm (ZLE) and new HRV indices (IDX PNA and IDX SNA ). ZLE clearly decomposed x RAW into x RSA and x R-f , and IDX PNA and IDX SNA identified parasympathetic activation and sympathetic inhibition during meditation, respectively. This study may be the first report on the identification of cardiac sympathetic inhibition during meditation (or slow breathing). Although the overlapping issue of the ZLE still remains, the ZLE was more robust than previous tachogram decomposition algorithms (Gauss, ARMAx and OSP) when RR dynamically fluctuates. Since the proposed approach requires only a single-lead ECG, we expect that it will be used in various fields, such as the internet of Medical Things (IoMT), digital healthcare, digital meditation, and digital therapeutics.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request. Table 3. Mean ± standard deviation of HRV indices (HF, IDX PNA , LF/HF, SB U , and IDX SNA ). Bold text represents a statistically significant difference between neutral and other affective states (p < 0.000005). HF and IDX PNA are the HRV indices related to vagal tone, and other indices (LF/HF, SB U and IDX SNA ) indicate sympathetic tone. In meditation state, IDX PNA s of ARMAx and ZLE showed significant increases, and statistically significant decreases were observed in SB U s and IDX SNA s of ZLE.